# Extreme weather events do not increase political parties' environmental attention

#'   This scirpt
#'    > performs PanelMatch models
#'    > creates one graphs for green party representation robustness


# Tim Wappenhans, António Valentim, Heike Klüver, and Lukas F. Stoetzer
# First: 18-02-2024
# Last: 23-02-2024



# PACKAGES --------------------------------------------------------------------
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  tidyverse,
  PanelMatch,
  panelView,
  ggpubr,
  ggrepel,
  cowplot,
  here
)


# DATA WRANGLING FOR PANELMATCH ------------------------------------------------
df.panelmatch <- read_csv(here("data/data_output/press_weekly.csv"))  |> 
  mutate(
    
    # only integers allowed 
    parlgov_id = as.integer(parlgov_id),
    t = as.integer(t),
    country = unclass(as.factor(country_name)),
    family = unclass(as.factor(family_name_short))
    
  ) |> 
  
  select(
    
    # Covs
    parlgov_id, t, family, country, cabinet_party,
    
    # Treatments
    treat, treat_flood, treat_storm, treat_fire, treat_temp,
    
    # Outcome
    issue_7, pr_total
    
  ) |> 
  
  mutate(
    
    outcome = issue_7 / pr_total,
    outcome = ifelse(is.nan(outcome), 0, outcome) 
    
  ) |> 
  
  as.data.frame() |> 
  
  filter(!is.na(parlgov_id)) 



# Check if there's a green party in each country*week
df.greens_in_parl <- df.panelmatch |> 
  mutate(green = ifelse(family == 5, 1, 0)) |> 
  group_by(country, t) |> 
  summarize(greens_in_parl = mean(green)) |> 
  mutate(greens_in_parl = ceiling(greens_in_parl))
  

# join indicator back to df.panelmatch
df.panelmatch <- left_join(df.panelmatch, df.greens_in_parl,
                           by = c("country", "t"))

# mutate outcomes

df.panelmatch <- df.panelmatch |> 
  mutate(
    attention_greens = ifelse(greens_in_parl == 1, outcome, NA),
    attention_no_greens = ifelse(greens_in_parl == 0, outcome, NA)
    )


# Globally define leads and lags, iterations
lag <- 6
lead <- 6
n_iterations <- 1000


# create empty tibble for data viz
gg.did <- data.frame(
  coef = numeric(), 
  se = numeric(),
  model = numeric()
)

# LOOPS TO CREATE ONE GRAPH ACROSS ALL MODELS ----------------------------------
models <- c("w/ Greens",
            "w/o Greens")
for (model in models) {
  if (model == "w/ Greens") {
    df.tmp <- df.panelmatch |> 
      mutate(outcome = attention_greens)
  } else if (model == "w/o Greens") {
    df.tmp <- df.panelmatch |> 
      mutate(outcome = attention_no_greens)
  } 
  
  PM.results <- PanelMatch(lag = lag, 
                           lead = 0:lead, 
                           data = df.tmp,
                           time.id = "t", 
                           unit.id = "parlgov_id",
                           treatment = "treat", 
                           outcome.var = "outcome",
                           match.missing = TRUE,
                           refinement.method = "mahalanobis",
                           covs.formula = ~ I(lag(outcome, 1:lead)),
                           exact.match.variables = "family",
                           size.match = 10, 
                           qoi = "att",
                           forbid.treatment.reversal = FALSE)
  
  
  ## Treatment
  PE.treat <- PanelEstimate(sets = PM.results, 
                            data = df.tmp,
                            se.method = "bootstrap",
                            number.iterations = n_iterations)
  
  # Store Coefs and SEs
  gg.treat <- summary(PE.treat, 
                      verbose = F, 
                      bias.corrected = F) |> 
    as.data.frame() |> 
    mutate(model = model,
           period = row_number() - 1) |> 
    rename(low_ci = `2.5%`,
           up_ci = `97.5%`) |> 
    select(-std.error)
  
  ## Placebo
  PE.placebo <- placebo_test(PM.results, 
                             data = df.tmp, 
                             number.iterations = n_iterations)
  
  # placebo estimates  
  placebo_estimates <- PE.placebo$estimates |> as_tibble() |> 
    mutate(period = row_number() - (lead + 1)) |> 
    rename(estimate = value) 
  
  
  # placebo CIs, percentiles from bootstrapping
  percentiles <- PE.placebo$bootstrapped.estimates |> as_tibble() |> 
    summarize(across(everything(), ~ quantile(., probs = c(0.025, 0.975)))) 
  
  percentiles_longer <- percentiles |> 
    mutate(quantile = row_number(),
           quantile = case_when(quantile == 1 ~ "low_ci",
                                quantile == 2 ~ "up_ci")) |> 
    pivot_longer(cols = starts_with("t-"), names_to = "period", values_to = "value") %>%
    pivot_wider(names_from = "quantile", values_from = "value") |> 
    mutate(period = substr(period,2,3),
           period = as.numeric(period))
  
  # join placebo estimates and placebo CIs
  gg.placebo <- left_join(placebo_estimates, percentiles_longer,
                          by = "period") |> 
    mutate(model = model)
  
  # bind placebo to treatment
  gg.tmp <- rbind(gg.placebo, gg.treat)
  
  # Add to main viz DF
  gg.did <- rbind(gg.did, gg.tmp)
  
}



# PLOTS ------------------------------------------------------------------------
# globally define settings
pre_color <- "#4885C1"
post_color <- "#AE3A4E"

pre_fill <- "#4885C1"
post_fill <- "#AE3A4E"


pre_shape <- 17
post_shape <- 16

# create indicator
gg.did <- gg.did |> 
  mutate(pre = ifelse(period < 0, "pre", "post"))



# data wrangling for plot
gg.did <- gg.did |> 
  mutate(
    
    # make it percentage points
    estimate = estimate * 100, 
    low_ci = low_ci * 100, 
    up_ci = up_ci * 100

    )


gg.did |> 
    ggplot() +
    
    # Gray pre treatment shade
    geom_rect(aes(xmin = -Inf, xmax = -0.5, ymin = -Inf, ymax = Inf),
              fill = "#e0e0e0", alpha = 0.2) +
    
    # CIs
    geom_linerange(aes(x = period, ymin = low_ci, ymax = up_ci, color = pre),
                   alpha = .7,
                   linewidth = 0.75) + 
    
    # Coef
    geom_point(aes(x = period, y = estimate, shape = pre, color = pre),
               size = 3) + 
    
    # fixed lines
    geom_hline(yintercept = 0, linetype = 'dotted') + 
    
    # look
    ylab("Percentage change in attention") +
    xlab("Weeks relative to extreme weather event") + 
    scale_x_continuous(breaks = gg.did$period, labels = gg.did$period) +
    facet_wrap(~ model, 
               scales = "free", 
               ncol = 2) +
    scale_color_manual(values = c("pre" = pre_color,
                                  "post" = post_color)) +
    scale_shape_manual(values = c("pre" = pre_shape, 
                                  "post" = post_shape)) +
    scale_fill_manual(values = c("pre" = pre_fill, 
                                 "post" = post_fill)) +
    guides(shape = "none",
           fill = "none",
           color = "none") +
    
    coord_cartesian(xlim = c(min(gg.did$period), max(gg.did$period))) + 
    
    # theme
    theme_bw() +
    theme(panel.grid = element_blank(),
          plot.title = element_text(hjust = 0.5)) +
  facet_grid(~ model)

# EXPORT -----------------------------------------------------------------------
ggsave(here("results/graphs/fig_ex_data_04.pdf"),
       width = 7,
       height = 5)